X-ray absorption Debye- Waller factors from ab initio molecular dynamics 
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An ab initio equation of motion method is introduced to calculate the temperature-dependent 
mean square vibrational amplitudes which appear in the Debye- Waller factors in x-ray absorption, 
x-ray scattering, and related spectra. The approach avoids explicit calculations of phonon-modes, 
and is based instead on calculations of the displacement-displacement time correlation function 
from ab initio density functional theory molecular dynamics simulations. The method also yields 
the vibrational density of states and thermal quantities such as the lattice free energy. Illustrations of 
the method are presented for a number of systems and compared with other methods and experiment. 



I. INTRODUCTION 

Thermal vibrations give rise to exponentially damped 
Debye-Waller (DW) factors exp[— W(T)] in x-ray ab- 
sorption spectra (XAS), x-ray diffraction (XRD), and 
related spectra. For example, in x-ray absorption fine 
structure spectra (XAFS) W(T) w 2k 2 a 2 R {T) where 
a R = ([(w-R — uo) ■ R] 2 ) refers to the mean square relative 
displacement (MSRD) of a given bond R 1 k is the photo- 
electron wave number, and T the absolute temperature. 1 
In XRD, and similarly in neutron diffraction (ND) and 
the Mossbauer effect, a 2 (T) = ((u ■ k) 2 ) is the mean- 
square displacement of an atom along the momentum 
transfer vector k, u being the instantaneous displacement 
vector. Due to their strong variation with temperature, 
energy, and the geometrical structure of a material, ac- 
curate DW factors are crucial to a quantitative analysis 
of XAS; conversely, the lack of precise Debye-Waller fac- 
tors is one of the main limitations to accurate structure 
determinations from experiment, especially for coordina- 
tion numbers. 

Various methods have been developed for obtaining 
these DW factors. Phenomenological models, e.g., cor- 
related Einstein and Debye models 2 ' 3 , are widely used 
in fitting but are often only semi-quantitative. More 
generally, they can be calculated in terms of Debye- 
integrals over appropriate projected vibrational densi- 
ties of states (VDOS). In small molecules, explicit sums 
over modes can be used to calculate the VDOS<^£ Such 
sums can also be used for periodic solids, both for crys- 
tallographic Debye-Waller factors and other thermody- 
namic quantities^— For complex materials, however, cal- 
culating and summing over modes can be a computa- 
tional bottleneck. As an alternative, a Lanczos algorithm 
can be used to evaluate the VDOS, starting from a dy- 
namical matrix (or Hessian), that can be obtained ei- 
ther from force-field models, 9,1 " or first principles DFT 
calculations »ii At high temperatures, brute force classi- 
cal MD (or DFT/MD) methods can also be used to ob- 
tain moments of vibrational distribution functions^ but 
such methods can fail at low temperatures when quantum 
statistics dominate. First principles DFT methods can 
also be computationally demanding, especially in com- 
plex systems. 



In an effort to speed up the calculations, we present 
here a first principles approach based on an ab initio 
equation of motion (AEM) approach using DFT molecu- 
lar dynamics calculations of displacement-displacement 
time-correlation functions. The method is a general- 
ization of the equation of motion method 1 - for calcula- 
tions of the VDOS, which was adapted for calculations of 
Debye-Waller factors based on force- field models. 14 How- 
ever since accurate force-field models are not generally 
available, especially for complex molecules and solids, 
DFT or other ab initio methods are needed. Because 
o\ depends primarily on the vibrational structure in the 
local environment around a given bond R, the calcula- 
tions can be carried out using relatively small clusters of 
atoms, without the use of periodic boundary conditions 
or other symmetry considerations. Thus the approach is 
applicable to general aperiodic materials. 



II. EQUATION-OF-MOTION METHOD 
A. Formalism 

The theory used in the present study is a first princi- 
ples extension of the equation of motion approac h 13 i for 
calculations of the VDOS and thermodynamic quantities 
that can be expressed as Debye-integrals over the VDOS. 
Our ab initio equation of motion (AEM) extension builds 
in dynamical structure in terms of first principles DFT 
calculations for a general structure, but does not rely 
on explicit calculations of the dynamical matrix (DM). 
The technique builds in Bose-Einstein statistics, and al- 
lows one to calculate the DW factors and related thermal 
properties either in real-time or the frequency domain. 
The AEM method has a number of computational ad- 
vantages. It can be efficient even for large systems, since 
the method is local and diagonalization of huge matri- 
ces is not necessary. Also the computational time scales 
linearly with the size of a cluster. Anharmonic effects 
such as lattice expansion can be added using a cumulant 
expansion. 11 

Our AEM method is based on calculations of the 
displacement-displacement correlation function in real 
time, using solutions of the 3N coupled Newton's equa- 
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tions of motion with DFT/MD methods. Such corre- 
lation functions are Fourier transforms of projected vi- 
brational densities of states (VDOS), which are defined 
uniquely by the initial conditions. Physically the VDOS 
can be interpreted as the "sound" of a lattice "plucked" 
along a given set of initial displacements. Here N is the 
number of atoms in the system which is centered within 
the region of interest and typically a few near-neighbors 
in radius. Regarding the total lattice potential energy $ 
of the crystal lattice as a function of the local atomic dis- 
placements Ui from their thermal-equilibrium positions 
Ri, and making use of a quasi-harmonic approximation, 
the equations of motion can be written as 

^2 = - J! D ia,kf3Qkl3, (1) 

fc/3 

with given initial displacements Wj(0), (i = 1 . . . N), and 
zero initial velocities Uj(0) = 0. Here Qi = ui \JMi denote 
reduced displacements at site i where Mi is the atomic 
mass, and D ia ^p = ^i a ,kpl 'VMiMk is the dynamical 
matrix of order 3N x 3N. The matrix <f>i a k8 consists of 
second derivatives of the potential energy with respect 
to the atomic displacements Ui a and Uup, where i, k are 
atomic sites and a, ft = (x,y,z). Formally, the reduced 
displacement vectors Qi can be expanded in normal co- 
ordinates q\ and eigenmodes A as 

Q 4 =]Te*(A) gA . (2) 

A 

Substituting this relation into Eq. ([1} , leads to a standard 
eigenvalue problem for the normal modes, 

£,a (A) = ^ D ia,k/3 £fe/3 (A) . (3) 

After evaluating the thermal averages using Bose- 
Einstcin statistics, one obtains for the normal coordinates 

"l(qx) 2 = (n(w A ) + = coth (4) 

In applying these results for calculations of interest 
here, it is convenient to define a normalized displacement 
state 

\Q(t)) = \Q 1 (t),Q 2 (t),...Q N (t)). (5) 

For example, for the MSRD for a given near-neighbor 
bond (0, -R)r^ the initial displacement state \Qr(0)} has 
Qo(0) =-vW^i Qr(°) = +y/liR/M R R, and oth- 
erwise Qi = 0, where (ir — (1/Mr + 1/Mq)" 1 is the 
reduced mass. A frequency domain expression for the 
MSRD can then be obtained from Eqs. (J2JH]) and sum- 
ming over all modes, i.e., 

a% = ([(u R - fib) • R?) 

= J^ T \(MQn(Q))\\ oth ^ (6) 

2fi R ^ lj x 2 

= 7) — / — PR\u) coth-—. (7) 
z Hr. Jo u z 



Here 

p«H^|(A|Q fl (0))| 2 ^-^A) (8) 

A 

is the projected VDOS contributing to relative vibra- 
tional motion along R and /3 = 1/fcgT. The maximum 
frequency oj max in Eq. j7]) can be estimated from the 
relation Lo max > z^Jk/ fiR where z is the coordination 
number and k is the near-neighbor force constant. 

In order to obtain an equivalent time-domain expres- 
sion for the VDOS, we calculate the cosine-transform of 
the displacement-displacement time-correlation function 
(Q R(t)\Q r(0)) with an ad hoc exponential damping fac- 
tor that limits the maximum time t max of the integration 
and MD runs, 

Pr (oj) = - / (Q R (t)\Q R (0)) cos ute~ Et2 dt (9) 
= ^|(A|Qh(0))| 2 <5a(^-^a). (10) 

A 

Thus as a consequence of the damping factor, the pro- 
jected VDOS of Eq. (TIT))) is broadened by narrow <5-like 
functions of width A typically chosen to be about 5% of 
the bandwidth. This broadening also smooths the oth- 
erwise discrete spectrum of the finite system used, but 
has practically no effect on integrated quantities. The 
spectral width A is determined by the cutoff parameters 
£ = 3/ t2 nax an d tmax = \f§/{u m axA). These cutoff pa- 
rameters also focus on the local environment by cutting 
off long distance behavior. The time-correlation function 
in Eq. © is 

(Qfl(i)IQfl(O)) = Q^ a {t)Q la {Q). (11) 

where ur is the number of non- vanishing displacements 
m \Qr)- Instead of using the 2nd order differential equa- 
tions in Eq. ([]}, in our approach the displacement state 
vector \Qii(t)) is determined by integrating the equations 
of motion numerically using Velocity- Verle t 15 ' 16 molecu- 
lar dynamics with initial conditions as in \Qr(0)}, 

Ri(t + At) = Ri(t) +Vi(t)At + i a^At 2 , (12) 
Vi(t + At) = Vi(t) + - [ai{t)+ai(t + At)]At,(13) 

where Ri(t), Vi(t) and di(t) are, respectively, the instan- 
taneous position, velocity and acceleration of atom i. The 
acceleration Oj(i) = fi(t)/Mi and fi(t) is the force on 
atom i. The Hellmann-Feynman theorem ensures that 
the forces can be calculated as the expectation value of 
the analytical derivative of the Hamiltonian with respect 
to the nuclear positions. This algorithm is efficient since 
an explicit calculation of the dynamical matrix at each 
time-step is not necessary. 
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Finally a real time expression for the MSRD can be 
obtained by substituting Eq. © for Pr(uj) into ([7]) and 
evaluating the Fourier transform, yielding 



sharp resonances, due to the presence of poles in Eq. (|15l) . 
The coefficients can be obtained by solving the system 
of linear equations 



4.cn 



In 



dt(Q R (t)\Q R (0)) 



(2smh-)- 



(14) 



Eqs. (JT]), © and (fT4)l are the key formulas used in 
our AEM calculations. Throughout this work results 
obtained with Eqs. © and © will be labeled AEM- 
FT, while those obtained with Eq. (fl4|) will be labeled 
AEM-RT. The form of Eq. flU) shows that, it is not 
essential to determine the VDOS pr(uj) as an interme- 
diate step, and hence that cr R (T) can be calculated di- 
rectly from the corresponding displacement-displacement 
autocorrelation function. Note that in the time do- 
main the analog of the Bose-Einstein weight factor is 
-ln[2sinh7rt/^?i] = nt/(3h - In [exp (2wt/ph) - 1]. At 
long times when 7i/27ri < fceT the weight factor is nega- 
tive and reduces to — Ttt/j3h at high temperatures and to 
— In (2wt/{3h) at low. Due to the exponential damping, 
the net time integration limit t max is usually several vi- 
brational cycles and typically requires about 25-35 time- 
steps per cycle for accuracy to a few percent. In addition, 
the singular behavior of the integrands in Eq. © and 
(fT4"|) must be handled with care. This is especially im- 
portant at low temperatures due to zero-point motion. 
Thus in the time-domain, we further stabilize the long 
time behavior by convolving the time-correlation func- 
tion with the inverse Fourier transform of a smoothed, 
low- frequency cutoff function Q(uj — w c ), where uj c is an 
appropriate cutoff frequency. In the frequency domain 
in Eq. © we replace the very low-frequency region with 
a similar cutoff or a Debye-model chosen to fit the very 
low frequency behavior of pr(u>). All the integrals in our 
implementation of the AEM method are evaluated us- 
ing the trapezoidal rule, which is appropriate for highly 
oscillatory integrands. 



B. 



Maximum Entropy Method 



Since the MSRDs are obtained from Debye-integrals 
over the VDOS, a precise determination of the spectra 
is not important, as long as the leading moments are 
accurate. Thus, as an alternative approach, the projected 
density of states can be obtained approximately using the 
Maximum Entropy Method (MEM)i^ In this approach 
the VDOS is approximated as 



M 

ate 

fc=i 



(15) 



where At is the sampling interval in the time domain and 
M is the desired order of the approximation. The MEM 
approach is well suited to represent phonon densities with 



-$^|3-A|°i = & (k = l,...M), (16) 



where 



N — l 



N 



^ ^(Qr{U)\Qr(P))(Qr(J^.i)\Qr(P))- ( 17 ) 



and N is the number of MD evolution steps. Although 
the MEM method can be less efficient than the direct FT 
approach, we find that it can be more stable in the low 
frequency region, since it is less sensitive to non-periodic 
trends in the time evolution. The reduced efficiency 
arises from the high order of approximation (M > 200) 
needed to achieve an accurate representation of pr(uj) at 
all frequencies. Throughout this work results obtained 
with Eqs. © and (T5]) will be labeled AEM-MEM. 



C. Multiple scattering aij 

The above real-time AEM method can also be used 
to calculate the MSRD <r| for a given XAFS multiple- 
scattering path j with rij legs. This MSRD corresponds 
to the mean-square fluctuation in the effective MS path 
length 



oj = ((SR) ) = 



(18) 



Here Aj = (Ru- + Ru+)/2, where Ru± represent the 
directional unit vectors between the site i and the sites 
i — 1 before and i + 1 after, along the multiple-scattering 
path j. In analogy with the single scattering results, we 
obtain expressions similar to Eq. © for cr| and Eq. © 
for Pj (w), but with the weights in mode A given by 



l<A|Qi(0))| a = 




(19) 



These weights can be interpreted as the normalized prob- 
ability that an initial displacement state \Qj(0)}, corre- 
sponding to a multiple-scattering path stretch, is in vi- 
brational mode | A). Thus the initial displacements in the 



state \Qj(0)) are Q 
the inverse reduced mass is 



Here 



1 " J 1 

a Mi 



I A, 



(20) 



i=l 



which is defined so that (Q j (0)\Q j (0)) — 1 and Pj(to) is 
normalized. 
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D. Other Dynamical Properties 

Other dynamical properties can be obtained similarly, 
by generalizing the seed-state \Qr(0)) appropriately^ 
For example, when the seed state is defined as a single- 
atom displacement, the resulting correlation function 
yields the mean square atomic displacements u 2 (T) in 
x-ray scattering DW factors. Also, when all symmetry 
unique Cartesian atomic displacements are added, one 
obtains the total VDOS per site pr{w). This permits 
calculations of thermodynamic functions such as the vi- 
brational free energy per site^ 



F(T) = 3k B T / duj In 



2 sinh 



pr(u), (21) 



where fee is the Boltzmann constant. Finally, if the 
\Qr(Q)) seed state is initialized with atomic displace- 
ments perpendicular to R instead of parallel to it, 
we can generate the mean-square transverse displace- 
ment a^_(T), which provides a correction to the lattice 



expansion 



li 



E. Computational Details 

The micro-canonical (i.e., NVE) ensemble MD simula- 
tions for the applications presented here were done using 
VASr»i£ for the crystalline systems and siest a 19 ' 20 for the 
Zn-imidazole complex. These codes were chosen on the 
basis of efficiency, although in principle, any program ca- 
pable of NVE dynamics can be interfaced with the AEM 
codes used in this work. The VASP simulations used stan- 
dard ultrasoft pseudopotentials, and were optimized for 
efficiency in MD runs. The Ge calculations used a 2x2x2 
fc-point grid with a plane- wave cutoff of 105 eV, while for 
ZrW 2 8 the grid was 4x4x4 and the cutoff was 297 eV. 
The siesta calculations used Troullier-Martins norm- 
conserving pseudopotentials — and standard double- £ ba- 
sis sets with a single polarization function (DZP). The 
confinement-energy shift defining the numerical atomic 
orbitals was 10 meV. Finally, the Hartree and exchange- 
correlation potentials were represented on a real-space 
grid with a plane- wave-equivalent cutoff of 120 Ry within 
a (18.4 A) 3 cell. Both crystalline and molecular simu- 
lations used the PBE functional. 2 — We have previously 
shown that the choice of exchange-correlation functional 
plays an important role in obtaining accurate MSRDs 
for metallic systemsAi However, here we only focus 
on non-metallic and molecular systems, for which the 
PBE functional yields reasonable accuracy compared to 
experiment^ 



F. Efficiency Considerations 

The efficiency of the AEM method depends on three 
factors: 1) The number of individual MSRDs that need to 



be computed, 2) the minimum and maximum frequencies 
that contribute to the VDOS, and 3) the quality of the ab 
initio MD. First, if a large number of MSRDs is needed, 
the computation of the full DM may be preferable since it 
yields all necessary DW factors with minimal additional 
effort. However, in most XAFS analysis only a handful 
of local DW factors need to be known accurately while 
those for more distant shells can approximated roughly 
using correlated Debye or Einstein models. For example, 
in the case of the coordination shell around a metallic 
center in a complex biomolecule the AEM approach can 
provide an efficient alternative to the Lanczos DM ap- 
proach. Second, if a given MSRD has similar contribu- 
tions from low and high frequency modes, the MD must 
have a short enough time-step to accurately represent 
the high frequency (25-35 steps per cycle) and a total 
run time with sufficient cycles of the low frequency (4-8 
cycles). Third, the AEM approach can take advantage 
of efficient implementations of DFT energies and forces 
such as those used here, without relying on analytic sec- 
ond derivatives needed in the Lanczos DM approach or 
the equations of motion in Eq. (JT|). 

Of the applications presented here, results for Ge and 
Zn +2 -tctraimidazole can be more efficiently treated us- 
ing the Lanczos DM approach. In the case of Ge this is 
due to the simplicity of the unit cell. In the case of Zn +2 - 
tetraimidazole, first there are a relatively small number of 
modes and second the modes cover a broad range of fre- 
quencies that would require small time-steps and a long 
total simulation time to represent accurately. On the 
other hand, the zirconium tungstate (ZrW20s) system, 
illustrates the definite advantage of the AEM approach 
for complex systems, since only a handful of MSRDs are 
needed for XAFS, while the unit cell contains hundreds 
of atoms. Based on our experience with the DM Lanczos 
approach, we estimate that the AEM approach would be 
nearly two orders of magnitude faster than a dynamical 
matrix calculation. 



III. APPLICATIONS 
A. Germanium 

As a relatively simple test case, the AEM was ap- 
plied to a crystalline germanium system using an 64-atom 
supercell generated by repeating 2x2x2 times the dia- 
mond cubic cell, with the experimental lattice constant 
of 5.6575 A. The MD simulations used a 2 fs timestep 
and a total simulation time of 4.5 ps. The initial struc- 
ture was generated by introducing a 4.8% bond stretch 
to one of the nearest neighbor pairs in the cell. 

The correlation function resulting from the Velocity- 
Verlet time evolution is shown in Fig. [TJ As expected, 
the oscillations are dominated by a single mode with a 
period of about 117 fs, associated with the Ge-Ge opti- 
cal mode stretch. This dominant behavior can also be 
observed in the VDOS shown in Fig. [21 where the opti- 
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FIG. 1: (Color online) Displacement-displacement correlation 
function for the nearest neighbor Ge-Ge bond, with and with- 
out a damping factor e = 7 x 1CP 7 fs~ 2 , obtained from a 
constant energy molecular dynamics simulation. 
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FIG. 2: (Color online) Phonon density of states projected on 
the nearest neighbor Ge-Ge interaction calculated with the 
AEM-MEM and AEM-FT approaches, and for comparison, 
the broadened Lanczos DM results. 



cal modes are centered at about 8.5 THz. An integra- 
tion time of about 2 ps is adequate to obtain phonon- 
spectra with a spectral broadening of about 5%. The 
centroid of the VDOS is located at about 8 THz, in good 
agreement with Einstein models for the nearest neighbor 
single-scattering path with an Einstein frequency of 7.55 
THz. 23 It should be noted that although the integration 
time for optical mode is well above that needed for con- 
vergence, the net integration time for lower frequencies 
around 5 THz is just adequate. Due to the singular be- 
havior in Eq. (j6|), an adequate time integration for the 
lower frequency components is essential, and is especially 
important at low temperatures for some of the systems 
discussed in the next section. 

The MSRDs calculated for the nearest neighbor Ge-Ge 
bond are shown in Fig. [3J The agreement with experi- 
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FIG. 3: (Color online) Nearest neighbor Ge-Ge MRSD 
calculated with the AEM-FT, AEM-MEM and AEM- 
RT approaches, and for comparison, Lanczos DM and 
experimental results. The experimental results are shifted 
as in Ref . [Til. 



ment is quite good, with an average error of 4% for the 
AEM-FT approach and 2% for the AEM-MEM approach. 
For comparison, the DM-Lanczos approach has an aver- 
age error of 2%. Fig. [3] also shows the results obtained 
with the real-time approach of Eq. (|14[) and a frequency 
cutoff of 1.7 THz as in the FT and MEM approaches. As 
expected, given the formal equivalence between Eq. (fT4"|) 
the FT approach with an intermediate calculation of the 
VDOS, the results are nearly identical. 

To explore the accuracy and efficiency of the AEM-FT 
and AEM-MEM approaches, we have also integrated the 
correlation function both for shorter times and for larger 
time steps. For Ge we find that the total integration 
time can be reduced to about 1 ps without significant 
loss of accuracy. This corresponds to about 10 periods of 
the 8.5 THz dominant frequency. For integration times 
of about 500 fs the mean error for the MSRD increases 
to 8% for the FT approach and to 16% for MEM. From 
the point of view of the length of the time step, both 
the FT and MEM approaches are extremely resilient. In 
both cases the mean errors for the Ge MSRD remain 
constant with time steps up to 24 fs. This corresponds 
to approximately five samples per period of the 8.5 THz 
frequency. Such large time steps, however, might not be 
feasible within the MD simulation itself due to loss of 
energy conservation in the Verlet algorithm. 

As an example of other dynamical quantities that can 
be obtained with the AEM approach, Fig. @] shows the 
total phonon density of states for Ge calculated with the 
FT and MEM approaches. For comparison broadened 
dynamical matrix Lanczos and experimental^ results are 
also included. This VDOS was obtained by applying a 
single atomic displacement along the x axis, as described 
in III CI and by propagating as for a 2 for 4.5 ps. Over- 
all, the centroid of the DOS is accurately reproduced by 
all methods: The centroid of the experimental DOS is 
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FIG. 4: (Color online) Total phonon density of states 
for Ge calculated with the AEM-MEM and AEM-FT ap- 
proaches, and compared to results from the Lanczos DM and 
experiment.— 

located at 5.8 THz, while the FT and MEM approaches 
place it at 6.0 and 5.7 THz, respectively. The spread (i.e., 
2nd moment) of the DOS is also well reproduced with 
the FT and MEM, giving 2.8 and 2.9 THz, respectively, 
versus 2.6 THz in the experiment. Finally, all methods 
reproduce the positions and weights of main features of 
the experimental VDOS quantitatively. On average the 
positions of the peaks deviate by at most 0.4 THz (i.e., 
about 4% of full bandwidth) and the relative weights are 
within 5% of those observed in experiment. 

The accuracy of the total VDOS can also be gauged 
by comparing with the experimentally measured atomic 
MSD u 1 for Ge. Fig. [3] shows the MSD computed using 
the total VDOS shown in Fig. H The AEM results are 
in excellent agreement with those obtained with the full 
DM Lanczos approach and in good agreement with the 
available experimental results^ except at low tempera- 
tures. 



B. Zn -tetraimidazole 

As an example of a complex molecule, Zn +2 - 
tetraimidazole was simulated using the full structure 
shown in Fig. [51 This structure was optimized in siesta 
and one of the equivalent Zn-N bonds was distorted with 
a 3.4% bond stretch. The MD simulations used a 3 fs 
timestep and a total simulation time of 3.9 ps. Given 
its large number of degrees of freedom, the dynamics of 
Zn +2 -tetraimidazole are significantly more complicated 
than those of Ge. This can be seen in the correlation 
function shown in Fig. which exhibits a superposition 
of several modes. The dominant contributions can be an- 
alyzed by examining the VDOS in Fig. [8] The DM ap- 
proach exhibits three dominant frequencies at 5, 13 and 
25 THz, which contribute 32, 18 and 24%, respectively, 
of the MRSD value. It is interesting to note that the 
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FIG. 5: (Color online) Mean square atomic displacement 
for Ge calculated with the AEM-FT approach using a sin- 
gle atomic displacement, and for comparison, Lanczos DM 
and experimental results. 




FIG. 6: (Color online) Structure of Zn +2 -tetraimidazole. 



weight of the associated poles is 9, 13 and 31%, further 
highlighting the importance of the correct representation 
of the low frequency modes. In principle, the Zn-N path 
should be dominated by low frequency Zn-ligand tetrahe- 
dral modes. Loeffen et alA find that these modes appear 
at about 6.5 THz, in fair agreement with our principal 
contribution at 5 THz. Although the VDOS calculated 
with the AEM-FT and AEM-MEM approaches are in 
good agreement with each other, they have small differ- 
ences with respect to the Lanczos DM VDOS. For in- 
stance, the mode at 25 THz is blueshifted about 2 THz 
in the real-time approaches. Fig. [9] shows that the agree- 
ment between the MSRDs calculated from the different 
VDOS is quite good. At 8% error, the theoretical results 
are less accurate than those obtained for Ge. They are, 
however, still within the error margins of the available 
experimental value at 20K. The larger error is likely due 
to the quality of the basis set used in the SIESTA calcu- 
lations. 
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FIG. 7: (Color online) Displacement-displacement correla- 
tion function for the nearest neighbor Zn-N interaction in 
Zn +2 -tetraimidazole with and without a damping factor e = 
6 x 10 -7 fs -2 , obtained from a constant energy molecular 
dynamics simulation. 
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FIG. 8: (Color online) Phonon density of states projected on 
the nearest neighbor Zn-N interaction in Zn +2 -tetraimidazole 
calculated with the AEM-MEM and AEM-FT approaches, 
and for comparison, the broadened Lanczos DM results. 
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FIG. 9: (Color online) Nearest neighbor Zn-N MRSD in 
Zn +2 -tetraimidazole calculated with the AEM-MEM and 
AEM-FT approaches, and for comparison, Lanczos DM and 
experimental results. 




FIG. 10: (Color online) Structure of 2x2x2 supercell of zir- 
conium tungstate ZrW208 (Zr: light blue, W: dark blue, O: 
red). 



C. ZrW 2 O s 

Our final example is zirconium tungstate (ZrM^Os), a 
ceramic that exhibits negative thermal expansion (NTE) . 
This system is quite challenging, having a complex unit 
cell that puts the calculation of the DM for the Lanczos 
approach beyond the reach of our current implementa- 
tion and computational capabilities. Here we have ap- 
plied the AEM approach to a 352-atom supercell (Fig. 
mijl made of 2x2x2 repetitions of the unit cell. The sim- 
ulations used the experimental unit cell lattice constant 
of 9.1546 A and a timestep of 4 fs, for a total simula- 
tion time of 1.5 ps. ZrW20s has several interactions of 
interest, including Zr-Zr, W-W, W-0 and two inequiv- 
alent nearest-neighbor Zr-0 bonds with distances 2.03 



and 2.11 A. In principle, any of these interactions can be 
studied using the AEM approach. As a proof of principle 
here we study the a 2 of the shortest of the Zr-0 bonds by 
using an initial structure corresponding to a 3.8% bond 
stretch. 

For a Zr-0 distortion, the dynamics of Z1-W2O8 are not 
as complex as those observed for Zn +2 -tetraimidazolc. 
The correlation function (Fig. Hip is mostly dominated 
by a mode with a 40 fs period superposed on a mode 
with a period approximately three times longer. Visual 
inspection of the MD trajectory reveals that the 40 fs 
mode is associated principally with the longitudinal Zr- 
O stretch mode. These vibrational modes can be clearly 
seen at about 25 and 8 THz, respectively, in the VDOS 
shown in Fig. 1121 As in the previous examples, the agree- 
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FIG. 11: (Color online) Displacement-displacement correla- 
tion function for the short Zr-O bond nearest neighbor in- 
teraction, in ZrW2 08 with and without a damping factor 
e — 6x 10 -6 fs -2 , obtained from a constant energy molecular 
dynamics simulation. 
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FIG. 12: (Color online) Phonon density of states projected on 
the short Zr-0 bond nearest neighbor interaction of ZrW2 08, 
calculated with the AEM-FT and AEM-MEM approaches, 
and for comparison the experimental Raman spectrum.— 



merit between the AEM-MEM and AEM-FT approaches 
is very good. The agreement with the mode frequencies 
observed in the experimental Raman spectrum is also 
quite good. The FT and MEM VDOS show modes at 
approximately 7.7, 24.5, 27.7 and 31.7 THz, compared to 
the experimental peaks at 5.7-11.8, 23.8, 27.9 and 31.0 
THz. The 5.7-11.8 THz peaks are associated mostly with 
modes located on the tungstate ion and with some low 
frequency W0 4 modest The 23.8, 27.9 and 31.0 THz 
peaks correspond exclusively to asymmetric WO4 modes. 
It is interesting to note that the dynamics of this system 
are quite complex. Since the WO4 units are very stiff, 
the ZrC>6 units must rotate as the WO4 units translate.- 7 
Thus, a simple distortion of the Zr-0 bond is able to ac- 
tivate both the low and high frequency modes. 



FIG. 13: (Color online) MRSD of the shortest nearest neigh- 
bor Zr-O bond in ZrW 2 O s calculated with the AEM-FT, 
AEM-MEM and AEM-RT approaches, and for comparison 
experimental results. 



Fig. [13] shows the nearest-neighbor Zr-0 MSRD as a 
function of temperature. As expected, given the similar- 
ity of their VDOS, the FT and MEM values are in very 
good agreement. The direct integration RT approach 
also agrees well with the FT approach, at least for higher 
temperatures. The agreement with experimental is also 
quite good, with all theories falling within the experimen- 
tal error bars for most of the temperature range. The 
largest disagreement occurs in the 80-140K region where 
other MSRDs (W-W and Zr-Zr) are known to have an 
anomaly that is likely related to the NTEi22 



IV. CONCLUSIONS 

We have introduced an ab initio equation of mo- 
tion (AEM) method for calculations of the MSRDs a 2 , 
needed for Debye- Waller factors in x-ray absorption, x- 
ray scattering, and related spectra. The method is based 
on calculations of displacement-displacement time cor- 
relation functions from ab initio density functional the- 
ory molecular dynamics simulations, using the Velocity- 
Verlet time-evolution algorithm. Thus the approach 
avoids the need for explicit calculations of phonon-modes 
or the dynamical matrix. The AEM method builds in 
Bose-Einstein statistics and yields the vibrational den- 
sity of states (VDOS) as either cosine Fourier Trans- 
forms of displacement-displacement correlation functions 
or through the Maximum Entropy Method. The MSRDs 
and other thermal quantities such as the lattice free en- 
ergy, are obtained in terms of Debye-integrals over the 
VDOS. Alternatively, the MSRDs can be computed di- 
rectly from the correlation functions by using the time- 
domain counterpart of the Bose-Einstein weight factor. 
Application of the method to a number of systems show 
that the approach is computationally advantageous for 
large, complex systems, and is in quantitative agreement 
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with other methods and with experimental results. 
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